!> author: 左志华
!> date: 2022-07-16
!>
!> Load particle data from file <br>
!> 加载粒子信息
module sph_load_pif

    use, intrinsic :: iso_fortran_env, only: stderr => error_unit
    use seakeeping_string, only: to_string
    use seakeeping_filesystem, only: is_exist, operator(.join.)
    use sph_kinds, only: rk
    use h5part_m
    use sph_region, only: region_t, save_virtual_particles
    use sph_particle, only: particle_t
    use sph_save_pif, only: nstep
    use sph_math, only: const
    use sph_env, only: dim_ => dim, nums
    use sph_terminal, only: term

contains

    !> Load particle information from file <br>
    !> 装载粒子信息
    subroutine load_particle_information(dir, region, nums, go_on, go_on_step, hsml)
        character(*), intent(in) :: dir  !! 粒子信息文件路径
        type(region_t), intent(out) :: region  !! 粒子宏观变量域
        integer, intent(out) :: nums(5)  !! 粒子数目
        logical, intent(in) :: go_on  !! 是否继续
        integer, intent(in) :: go_on_step  !! 继续的步数
        real(rk), intent(out) :: hsml  !! 光滑长度

        if (go_on) nstep = int(go_on_step, 8)
        call region%load_pif_nml(dir.join.'pif.nml', nums, hsml)
        call h5part_load_pif(dir.join.'pif.h5part', region%particles, nums(4))
        region%particles%hsml = hsml
        allocate (region%virtual_particles%loc(dim_, nums(2)), &
                  region%virtual_particles%vel(dim_, nums(2)), &
                  region%virtual_particles%rho(nums(2)), &
                  region%virtual_particles%p(nums(2)))
        call save_virtual_particles(region%particles, region%virtual_particles, nums(2), nums(1))

    end subroutine load_particle_information

    !> 读取 hdf5 粒子信息
    subroutine h5part_load_pif(file, particle, ntotal)
        character(*), intent(in) :: file  !! 文件路径
        type(particle_t), intent(out) :: particle  !! 粒子宏观变量域
        integer, intent(in) :: ntotal  !! 粒子总数
        integer :: ierr, i
        integer(8) :: fid, status
        real(rk), allocatable :: z(:, :)

        fid = h5pt_openr(file)
        if (fid == 0) then
            call term%error('file may be in use by another program: '//file)
        end if

        allocate (particle%loc(dim_, ntotal), &
                  particle%vel(dim_, ntotal), &
                  particle%mass(ntotal), &
                  particle%rho(ntotal), &
                  particle%p(ntotal), &
                  particle%u(ntotal), &
                  particle%itype(ntotal), stat=ierr)
        if (ierr /= 0) then
            call term%error("failed to allocate space for a particle")
        end if

        allocate (particle%acc(dim_, ntotal), &
                  particle%drho(ntotal), &
                  particle%du(ntotal), &
                  particle%c(ntotal), &
                  source=0.0_rk, &     ! 初始化为 0
                  stat=ierr)
        if (ierr /= 0) then
            call term%error("failed to allocate space for a particle")
        end if

        ! ------------------------ 读取 0 时间步粒子数据 -------------------------- !
        status = h5pt_setstep(fid, nstep)

        allocate (z(ntotal, dim_), stat=ierr)
        if (ierr /= 0) then
            call term%error("failed to allocate space for a particle")
        end if

        status = h5pt_readdata(fid, "x", z(:, 1))
        status = h5pt_readdata(fid, "y", z(:, 2))
        do concurrent(i=1:nums(4))
            particle%loc(:, i) = z(i, :)
        end do

        status = h5pt_readdata(fid, "Vx", z(:, 1))
        status = h5pt_readdata(fid, "Vy", z(:, 2))
        do concurrent(i=1:nums(4))
            particle%vel(:, i) = z(i, :)
        end do

        status = h5pt_readdata(fid, "Mass", particle%mass)
        status = h5pt_readdata(fid, "Density", particle%rho)
        status = h5pt_readdata(fid, "P", particle%p)
        status = h5pt_readdata(fid, "Itype", particle%itype)

        status = h5pt_close(fid)

        ! ------------------------------- 检查 --------------------------------- !
        if (any(particle%rho(:) < const(1))) then
            call term%error("density is negtive: "//file)
        end if

    end subroutine h5part_load_pif

end module sph_load_pif

!> author: 左志华
!> date: 2022-07-16
!>
!> Load TOML config file <br>
!> 读取 TOML 配置文件
module sph_load_toml

    use seakeeping_string, only: to_string
    use seakeeping_filesystem, only: is_exist, operator(.join.)
    use sph_kinds, only: rk
    use sph_project_type, only: project_t
    use sph_save_pif, only: save_acc
    use sph_eos, only: c0, eos_type, get_sound_speed, H_SWL, B
    use sph_artifical_force, only: art_heat, art_visc
    use sph_density_summation, only: density_summation_type
    use sph_internal_force, only: viscosity
    use sph_external_force, only: gravity
    use sph_linked_list_nnps, only: estimated_capacity
    use sph_smoothed_kernel_function, only: skf
    use sph_navier_stokes_solver, only: nsave, dt
    use sph_env, only: has_rigid

contains

    !> Parse sph.toml <br>
    !> 获取项目信息
    subroutine parse_project_toml(project, working_dir)
        type(project_t), intent(inout) :: project  !! 工程结构体
        character(len=*), intent(in) :: working_dir  !! 工作目录
        character(:), allocatable :: sph_toml
        character(128) :: name, skf_type
        integer :: itime_integration, go_on_step
        logical :: go_on, adaptive_step_size
        real(rk) :: cfl
        namelist /sph/ name, dt, itime_integration, skf_type, go_on, go_on_step, &
            density_summation_type, nsave, save_acc, viscosity, art_visc, gravity, &
            adaptive_step_size, eos_type, CFL, c0, H_SWL, B, estimated_capacity, has_rigid

        project%path = working_dir
        skf_type = 'cubic_spline'
        sph_toml = project%path.join."sph-config.nml"
        open (1, file=sph_toml, action='read', status='old')
        read (1, nml=sph)
        close (1)

        project%name = trim(name)
        project%itime_integration = itime_integration
        project%go_on = go_on
        project%go_on_step = go_on_step
        project%adaptive_step_size = adaptive_step_size
        project%CFL = CFL
        skf%type = trim(skf_type)

        if (abs(c0 - 0.0_rk) < 1.0e-6_rk) then
            call get_sound_speed(H_SWL, c0)
        end if

    end subroutine parse_project_toml

end module sph_load_toml
